Prioritizing Potential Aquaculture
With the U.N. predicting that the global human population will reach approximately 9.7 billion people by the year 2050, access to basic resources (particularly food and water) will continue to become increasing pressing issues (United Nations 2024). The global assessment of the state of food security and nutrition in 2022 highlights the severity of this issue at current population levels (~8.2 billion):
“It is estimated that between 691 and 783 million people in the world faced hunger in 2022. Considering the midrange (about 735 million), 122 million more people faced hunger in 2022 than in 2019, before the pandemic (Food and Agriculture Organization of the United Nations 2023).”
This highlights the increasingly urgent need to address the lack of global food security. As a result, marine aquaculture has become the center of research; with many notable characteristics that make it a sustainable alternative to meet the growing demands for protein sources (Hall et al. 2011; Nash et al. 2021). Gentry et al. mapped potential locations that would be suitable for marine aquaculture around the world, and they found that global seafood demand could be met using less than 0.015% of the global ocean area (Gentry et al. 2017).
They identified potential locations based on human and environmental constraints. This included excluding areas; that would negatively impact the growth of bivalves (low phytoplanktonic food availability) and finfish (low dissolve oxygen), were too expensive (too deep >200m) to anchor farms, were already allocated to other uses, or were in high-density shipping areas (Gentry et al. 2017).
Based on this analysis, we aim to identify locations along the west coast of the U.S. that would be optimal for the establishment of a marine aquaculture farm for both oysters and various finfish species based on two environmental parameters that impact optimal growth:
Research has shown that the following conditions are required for optimal oyster growth:
Potential finfish species were selected based on the current commercial fishes managed by the NOAA West Coast operation (NOAA Fisheries: West Coast 2024). Highly migratory species (such as sharks) were removed from consideration resulting in a list of 26 potential finfish species. We then used FishBase to compile a table of minimum temperature, maximum temperature, minimum depth, and maximum depth for each species (FishBase Consortium 2024). If a temperature range was not explicity listed, the preferred minimum temperature and preferred maximum temperature (estimated from a model) were recorded (Froese 2020).
The first goal of this assessment, is to determine locations within the US west coast Economic Exclusion Zones (EEZ) that would be suitable for oyster cultivation. Then a generalizable workflow was developed to assess the suitability of the 26 potential finfish species previously identified. Finally, the results of the analysis are provided.
First, load the necessary packages for this analysis:
#this clears out the environment
rm(list = ls())
#load necessary data packages
library(terra)
library(sf)
library(tidyverse)
library(tmap)
library(here)
library(janitor)
library(knitr)
library(stringr)
library(tibble)
library(kableExtra)
library(RColorBrewer)
library(htmltools)All functions used in this workflow are outlined and defined in this section, and grouped into three subcategories:
This function can be used to create a tibble with all of the geometric properties (resolution, extent, origin, and coordiante reference system) for a list of spatial objects with the class “SpatRaster” or “sf object” (both terra and sf compatible). This is useful for comparing the spatial properties of the spatial objects in a workflow to ensure compatibility before processing begins.
The input is:
The output is:
#write a function to create a tibble with all the spatial properties for a list of spatial objects
extract_spatial_properties <- function(spatial_file_list) {
#create an empty tibble to hold the output
spatial_properties <- tibble(name = character(),
resolution = character(),
extent = character(),
origin = character(),
crs = character())
#write a for loop to extract the properties into the tibble for review
for (file_name in names(spatial_file_list)) {
#extract the current spatial object
obj <- spatial_file_list[[file_name]]
#since we have sf objects and SpatRasters, we must handle them differently:
##if the object is a SpatRaster
if (inherits(obj, "SpatRaster")){
#if the object is a raster stack (has more than one layer):
if (nlyr(obj) > 1) {
##for each layer of the raster stack
for (layer_name in names(obj)) {
#extract the layer
layer <- obj[[layer_name]]
#extract the spatial properties for this layer:
##extract the resolution (cell size) as a string with 5 decimals
resolution <- paste(round(res(obj), 5), collapse = ", ")
##extract the extent (raster size) as a string with 2 decimals
extent <- paste(round(ext(obj), 2), collapse = ", ")
##extract the origin (0,0) as a string with 7 decimals
origin <- paste(round(origin(obj), 7), collapse = ", ")
##extract the crs from the current object
crs <- terra::crs(obj)
#extract the EPSG ID part from the CRS string to make the CRS easier to read
crs_id <- sub('.*ID\\["EPSG",(\\d+)\\].*', 'EPSG:\\1', as.character(crs))
#add the new information to the tibble for raster properties
spatial_properties <- spatial_properties %>%
add_row(name = layer_name,
resolution = resolution,
extent = extent,
origin = origin,
crs = crs_id)
}
#if the object is a single raster:
} else {
#extract the spatial properties for this layer:
##extract the resolution (cell size) as a string with 5 decimals
resolution <- paste(round(res(obj), 5), collapse = ", ")
##extract the extent (raster size) as a string with 2 decimals
extent <- paste(round(ext(obj), 2), collapse = ", ")
##extract the origin (0,0) as a string with 7 decimals
origin <- paste(round(origin(obj), 7), collapse = ", ")
##extract the crs from the current object
crs <- terra::crs(obj)
#extract the EPSG ID part from the CRS string to make the CRS easier to read
crs_id <- sub('.*ID\\["EPSG",(\\d+)\\].*', 'EPSG:\\1', as.character(crs))
#add the new information to the tibble for raster properties
spatial_properties <- spatial_properties %>%
add_row(name = file_name,
resolution = resolution,
extent = extent,
origin = origin,
crs = crs_id)
}
##if the object is a sf object
} else if (inherits(obj, "sf")) {
#extract the spatial properties for this layer:
##extract the EPSG ID part from the CRS string to make the CRS easier to read
crs_id <- gsub('.*ID\\[\\\"(EPSG)\\\",(\\d+)\\].*', '\\1:\\2', as.character(crs))
#add the new information to the tibble for raster properties
spatial_properties <- spatial_properties %>%
add_row(name = file_name,
resolution = NA,
extent = NA,
origin = NA,
crs = crs_id)
}
}
#return the tibble with the extracted spatial properties
return(spatial_properties)
}This function can be used to verify that all the spatial objects in a list (both terra and sf compatible) have a pre-defined coordinate reference system (crs). This is necessary to ensure before any map algebra can be performed.
The input is:
The output is:
#write a function to verify that all the spatial objects have the same crs
check_crs <- function(obj, target_crs) {
#if the object is a "SpatRaster"
if (inherits(obj, "SpatRaster")) {
#returns a logical (TRUE/FALSE) if the crs matches the target_crs (TRUE) or not (FALSE)
return(terra::crs(obj) == target_crs)
#if the object is an "sf object"
} else if (inherits(obj, "sf")) {
#returns a logical (TRUE/FALSE) if the crs matches the target_crs (TRUE) or not (FALSE)
return(st_crs(obj)$wkt == target_crs) ##need to use wkt since terra and sf print crs differently
}
}This function can be used to check that the geometric properties (resolution, extent, origin, and crs) of two SpatRaster objects (only compatible with terra) match. This function is a modified version of the extract_spatial_properties function.
The input is:
The output is:
#write a modified version of the extract_spatial_properties function to check the spatial properties
##this is only compatible with raster objects (not sf objects)
check_spatial_properties <- function(spatial_file_list) {
#create an empty tibble to hold the output
spatial_properties <- tibble(name = character(),
resolution_match = character(),
extent_match = character(),
origin_match = character(),
crs_match = character())
#we will use the first file in the list as the reference object for verification
ref_obj <- spatial_file_list[[1]]
#now identify the spatial properties of the first object
##identify the reference resolution
ref_resolution <- paste(res(ref_obj), collapse = ", ")
##identify the reference extent
ref_extent <- paste(ext(ref_obj), collapse = ", ")
##identify the reference origin
ref_origin <- paste(origin(ref_obj), collapse = ", ")
##identify the reference crs
ref_crs <- terra::crs(ref_obj)
#write a for loop to extract the properties into the tibble for review
for (file_name in names(spatial_file_list)) {
#extract the current spatial object
rast_obj <- spatial_file_list[[file_name]]
#if the object is a raster stack (has more than one layer):
if (nlyr(rast_obj) > 1) {
##for each layer of the raster stack
for (layer_name in names(rast_obj)) {
#extract the layer
layer <- rast_obj[[layer_name]]
#extract the spatial properties for this layer:
##extract the resolution (cell size) as a string
resolution <- paste(res(rast_obj), collapse = ", ")
##extract the extent (raster size) as a string
extent <- paste(ext(rast_obj), collapse = ", ")
##extract the origin (0,0) as a string
origin <- paste(origin(rast_obj), collapse = ", ")
##extract the crs from the current object
crs <- terra::crs(rast_obj)
#check if resolution, extent, origin, and CRS match the reference raster
resolution_match <- ifelse(resolution == ref_resolution, "MATCH", "DOES NOT MATCH")
extent_match <- ifelse(extent == ref_extent, "MATCH", "DOES NOT MATCH")
origin_match <- ifelse(origin == ref_origin, "MATCH", "DOES NOT MATCH")
crs_match <- ifelse(crs == ref_crs, "MATCH", "DOES NOT MATCH")
#add the new information to the tibble for raster properties
spatial_properties <- spatial_properties %>%
add_row(name = layer_name,
resolution_match = resolution_match,
extent_match = extent_match,
origin_match = origin_match,
crs_match = crs_match)
}
#if the object is a single raster:
} else {
#extract the spatial properties for this layer:
##extract the resolution (cell size) as a string
resolution <- paste(res(rast_obj), collapse = ", ")
##extract the extent (raster size) as a string
extent <- paste(ext(rast_obj), collapse = ", ")
##extract the origin (0,0) as a string
origin <- paste(origin(rast_obj), collapse = ", ")
##extract the crs from the current object
crs <- terra::crs(rast_obj)
#check if resolution, extent, origin, and CRS match the reference raster
resolution_match <- ifelse(resolution == ref_resolution, "MATCH", "DOES NOT MATCH")
extent_match <- ifelse(extent == ref_extent, "MATCH", "DOES NOT MATCH")
origin_match <- ifelse(origin == ref_origin, "MATCH", "DOES NOT MATCH")
crs_match <- ifelse(crs == ref_crs, "MATCH", "DOES NOT MATCH")
#add the new information to the tibble for raster properties
spatial_properties <- spatial_properties %>%
add_row(name = file_name,
resolution_match = resolution_match,
extent_match = extent_match,
origin_match = origin_match,
crs_match = crs_match)
}
}
#return the tibble with the extracted spatial properties verification results
return(spatial_properties)
}This function can be used to generate a map (Tennekes 2018) from a SpatRaster (only comaptible with terra rasters). In this workflow, this function is used to generate a map with the “optimal” and “not optimal” land classifications for each EEZ region for each potential finfish species.
The input is:
The output is:
#write a function to generate a map for the species raster object
generate_maps <- function(rast_list, map_title){
#create a base map to be used from the first raster in the list
combined_map <- tm_shape(rast_list[[1]]) +
tm_raster(palette = c("lightgrey", "tomato3"),
alpha = 0.5,
labels = c("Not Optimal", "Optimal"),
title = "Suitability") +
tm_layout(main.title = "Species Name",
main.title.position = "center",
main.title.size = 1.5,
main.title.fontfamily = "Times New Roman",
main.title.fontface = "bold",
legend.outside = TRUE)
#for every object in the raster list create a map by iterating through the remaining rasters in the list
if (length(rast_list) > 1) {
for (i in 2:length(rast_list)) {
#extract the raster from the list
raster_layer <- rast_list[[i]]
#save the raster layer name
layer_name <- names(rast_list)[i]
#add the new layer to the map
combined_map <- combined_map +
tm_shape(raster_layer) +
tm_raster(palette = c("lightgrey", "tomato3"),
alpha = 0.5,
legend.show = FALSE)
}
}
#save the map title based on the name_list generated with the species scientific names
map_title <- map_name_list[i]
#make the map name the species_name + _map
map_name <- paste0("./maps/", map_title, "_map.png")
#save the tmap
tmap_save(combined_map, filename = map_name)
#return the map
return(combined_map)
}This function can be used to reclassify the geometries in a SpatRaster (only compatible with terra) based on an optimal paramter range. It generates a new SpatRaster that has three classifications:
The input is:
The output is:
#write a function to reclassify the geometries in a raster based on the optimal parameter range
suitable_locations <- function(raster_obj, min_parameter, max_paramter) {
#duplicate the raster that will be reclassified
mask <- raster_obj
#create a reclassification matrix to reclassify the raster based on the optimal parameter range
rcl <- matrix(c(-Inf, min_parameter, 1, #below the minimum value ("below_range")
min_parameter, max_paramter, 2, #in the optimal range ("optimal")
max_paramter, Inf, 3), #above the maximum value ("above_range")
ncol = 3, byrow = TRUE)
#use the reclassification matrix to reclassify the duplicated raster with the rcl
mask <- classify(mask, rcl = rcl)
#assign the group to a classification (not optimal or optimal)
levels(mask) <- data.frame(id = c(1, 2, 3),
cats = c("below_range", "optimal", "above_range"))
#rename the mask to have the name of the raster_obj + _mask
mask_name <- paste0(names(raster_obj), "_mask")
names(mask) <- mask_name
#return reclassified raster object
return(mask)
}This function is used to generate a new raster object from two optimized SpatRaster (only compatible with terra). This function is formated to work with SpatRasters that were generated from the previously defined suitable_locations function. It generates a new SpatRaster that has two classifications based on the optimal parameter range of both SpatRasters (raster_1 and raster_2):
The input is:
The output is:
#write a function to generate a new raster object from two optimized rasters (this function is formatted to work with rasters generated from the suitable_locations locations defined above)
optimal_locations <- function(raster_1, raster_2){
#duplicate raster_1 to serve as the combined raster
combined <- raster_1
#first assign any locations with missing data from either raster as "not optimal" group 1
combined[raster_1 == 2 | raster_2 == 2] <- 1
#assign locations where BOTH parameters are "optimal" to group 2
combined[raster_1 == 2 & raster_2 == 2] <- 2
#assign locations that are not identified as optimal in both rasters as "not optimal" to group 1
combined[!(combined[] == 2)] <- 1
#assign the categories a name (optimal or not optimal)
levels(combined) <- data.frame(id = c(1, 2),
cats = c("not_optimal", "optimal"))
#return reclassified raster object
return(combined)
}This function is used to rasterize the individual geometries in a sf object based on the name of a column in the dataframe. In this workflow, this was used to rasterize the region geometries from the EEZ dataframe.
The input is:
The output is:
#write a function to rasterize individual geometries based on the column that has the corresponding name
rasterize_geometries <- function(df, template_raster, name_column) {
#create a list of names the name column
name_list <- unique(df[[name_column]])
#create an empty list to store the individual layers
rast_list <- list()
#create a duplicated empty raster from the template raster to be the raster stack
rast_stack <- rast(template_raster, nlyrs = 0)
#create an empty list to hold the bboxes
bbox_list <- list()
#write a for loop to extract the geometry for each id and save it as a raster in the list
for(unique_name in name_list) {
#filter the df to the geometry that corresponds to the current name
geometry <- df[df[[name_column]] == unique_name, ]
#create a duplicated empty raster from the template raster to be the raster layer
rast_layer <- template_raster
#rasterize the extracted geometry
geom_rast <- rasterize(geometry, rast_layer)
#rename the raster with the name
names(geom_rast) <- unique_name
#add the new raster to the list
rast_list[[unique_name]] <- geom_rast
#create a boundary box from the geometry
bbox <- st_bbox(geometry)
#rename the bbox with the unique_name + _bbox
bbox_name <- paste0(names(unique_name), "_bbox")
names(bbox) <- bbox_name
#add the new bbox to the list
bbox_list[[unique_name]] <- bbox
#add the new raster as a layer to the raster stack
rast_stack <- c(rast_stack, geom_rast)
}
#return a list with the individual layers, the raster stack, and the bbox_list
return(list(rast_list = rast_list,
rast_stack = rast_stack,
bbox_list = bbox_list))
}This function is used to create a mask of a SpatRaster (only compatible with terra) for each layer in a SpatRaster stack. It will then also generate a map for each of the masked layers generated. This function is a modified version of the rasterize_geometries function. In this workflow, this function is used to create a masked SpatRaster for each region in the EEZ with the input raster (raster_obj_input) being the SpatRaster generated for each potential species. The function also CONTAINS the generate_maps function to create the maps for each layer.
The input is:
The output is:
#write a function to create a list with a raster stack, rasters, and a map for each region
extract_regions <- function(raster_obj_input, rast_stack_input, map_name_list) {
#create an empty raster to store the output of the for loop
rast_list <- list()
#create a duplicate of the raster object to be the masked layer
rast_obj <- raster_obj_input
#create a duplicated empty raster from the raster object to be the raster stack
rast_stack <- rast(rast_stack_input, nlyrs = 0)
for (region_name in names(rast_stack_input)) {
#extract the region layer from the raster stack
region_layer <- rast_stack_input[[region_name]]
#use the region_layer to mask the rast_obj
region_layer_mask <- terra::mask(rast_obj, region_layer)
#rename the raster with the region_name
names(region_layer_mask) <- region_name
#add the new raster to the list
rast_list[[region_name]] <- region_layer_mask
#add the new raster as a layer to the raster stack
rast_stack <- c(rast_stack, region_layer_mask)
}
#save the map title based on the name_list generated with the species scientific names
map_title <- map_name_list[i]
#create a map from the layer
map <- generate_maps(rast_list, map_title)
#return a list with the individual layers, the raster stack, and the map
return(list(rast_list = rast_list,
rast_stack = rast_stack,
map = map))
}This function is used to create a tibble that ranks the regional SpatRasters (only compatible with terra) generated for each of the species using the previously defined functions. The regions are ranked based on the percent of the total area that is classified as optimal for cultivation. In this workflow, this function ranks the regions from 1 (best EEZ region for cultivation) to 5 (worst EEZ region for cultivation) based on a calculated optimal percent_cover.
The input is:
The output is:
#write a function to extract the create a tibble with each species and the rank for their optimal regions
rank_regions <- function(rast_stack){
#create a tibble to hold the output of the for loop
region_rank <- tibble(species = character(),
region = character(),
optimal_area = numeric(),
total_area = numeric(),
percent_cover = numeric(),
rank = numeric())
#for every layer in the rast_stack
for (layer_name in names(rast_stack)){
#extract the region layer from the raster stack
layer <- rast_stack[[layer_name]]
#calculate the area (km^2) for each cell
cell_area <- cellSize(layer, unit = "km")
#calculate the area of all the cells with the "not_optimal" (group 1) classification
not_optimal_area_km2 <- sum(values((layer == 1) * cell_area), na.rm = TRUE)
#calculate the area of all the cells with the "optimal" (group 2) classification
optimal_area_km2 <- sum(values((layer == 2) * cell_area), na.rm = TRUE)
#calculate the total area in the region
total_area_km2 <- optimal_area_km2 + not_optimal_area_km2
#calculate the percent cover of optimal area
percent_cover <- (optimal_area_km2 /total_area_km2) *100
#add the new information to the tibble
region_rank <- region_rank %>%
add_row(species = NA,
region = layer_name,
optimal_area = optimal_area_km2,
total_area = total_area_km2,
percent_cover = percent_cover,
rank = NA)
}
#rank the regions with the highest optimal percent cover being ranked 1
region_rank <- region_rank %>%
mutate(rank = ifelse(percent_cover == 0, NA,
rank(-percent_cover, ties.method = "max"))) #if the percent cover is zero assign the rank NA, and if there is a tie, assign it the max rank to be conservative
#return a table with the rankings for each region and each species
return(region_rank)
}This section sets up the filepaths and files for later analysis:
Then, define the file paths for the 1. sea surface temperature, 2. bathymetry, and 3. economic exclusion zone data:
#define the pathway to access the average annual sea surface temperature (sst) files
sst_files <- list.files(here::here("data"),
pattern = "average_annual_sst.*\\.tif$",
full.names = TRUE)
#define the pathway to access the depth file
depth_file <- here::here("data", "depth.tif")
#define the pathway to access the economic exclusion zone (eez) files
eez_file <- here::here("data", "wc_regions_clean.shp")
#define the pathway to access the csv with the potential finfish species data
species_file <- here::here("data", "potential_finfish_species.csv")Now that we have defined the file paths, we can read in the necessary spatial data:
#read in the sea surface temperature (sst) rasters
sst_rasters <- terra::rast(sst_files)
#write a for loop to rename each layer to have the format sst_year
for (i in seq_along(names(sst_rasters))) {
#extract the year for the corresponding file, by extracting the last four digits
year <- str_extract(basename(sst_files[i]), "\\d{4}") #extracting the last four digits
#format the name of the current layer to be sst_year and save it as the new name for the layer
names(sst_rasters)[i] <- paste("sst", year, sep = "_")
}#read in the depth raster
depth_raster <- terra::rast(depth_file)
#rename the layer
names(depth_raster) <- "depth"#read in the economic exclusion zone (eez) raster
eez <- sf::st_read(eez_file)
#create a vector with region_id's that correspond to the regions based on the geographic order
region_ids <- c("Washington" = 1,
"Oregon" = 2,
"Northern California" = 3,
"Central California" = 4,
"Southern California" = 5)
#create a region_id column in the eez dataset that assigns the regions a number from north (washington) to south (southern_california)
eez <- eez %>%
mutate(region_id = region_ids[rgn])
#create a version of the dataframe with the capitalized region names still
eez_caps <- eez %>%
clean_names() %>%
rename(region = rgn) %>%
select(region_id, region, area_km2, geometry)
#tidy up the data
eez <- eez %>%
clean_names() %>% #clean up the column names
mutate(region = str_replace_all(rgn, " ", "_") %>%
str_to_lower()) %>% #create a region column with the name in snake_case
select(region, area_km2, geometry) #only keep the necessary columnsNext, we need to examine the geometric properties (resolution, extent, origin, and coordinate reference system (crs)) for each of the spatial objects that we are working with. We will need to transform any spatial objects that have mismatched properties.
We will start by identifying the geometric properties of each of the spatial objects with the extract_spatial_properties function:
#create a summary list with all the files
spatial_file_list <- list(sst_rasters = sst_rasters,
depth_raster = depth_raster,
eez = eez)
#identify the spatial properties of all the spatial files in the list
spatial_properties <- extract_spatial_properties(spatial_file_list)
#generate a kable from the spatial properties extracted above
spatial_properties_kable <- spatial_properties %>%
kable("html",
col.names = c("Name", "Resolution", "Extent", "Origin", "CRS"),
caption = htmltools::tags$div(style = "text-align: center; font-size: 20px;",
htmltools::tags$strong("Geometric Properties Identification Results"))) %>%
kable_styling(full_width = FALSE, font_size = 14) %>%
row_spec(row = 0, bold = TRUE) %>%
kable_classic(html_font = "Times New Roman")
#print the kable
spatial_properties_kable| Name | Resolution | Extent | Origin | CRS |
|---|---|---|---|---|
| sst_2008 | 0.04166, 0.04166 | ext(-131.98, -114.99, 29.99, 49.99) | -7.6e-06, -3.8e-06 | EPSG:9122 |
| sst_2009 | 0.04166, 0.04166 | ext(-131.98, -114.99, 29.99, 49.99) | -7.6e-06, -3.8e-06 | EPSG:9122 |
| sst_2010 | 0.04166, 0.04166 | ext(-131.98, -114.99, 29.99, 49.99) | -7.6e-06, -3.8e-06 | EPSG:9122 |
| sst_2011 | 0.04166, 0.04166 | ext(-131.98, -114.99, 29.99, 49.99) | -7.6e-06, -3.8e-06 | EPSG:9122 |
| sst_2012 | 0.04166, 0.04166 | ext(-131.98, -114.99, 29.99, 49.99) | -7.6e-06, -3.8e-06 | EPSG:9122 |
| depth_raster | 0.00417, 0.00417 | ext(-132, -114, 29, 50) | 0, 0 | EPSG:4326 |
| eez | NA | NA | NA | EPSG:4326 |
It appears that the coordiante reference system for the bathymetry data (depth_raster) and the economic exclusion zone data (eez) does not match that of the sea surface temperature data (sst_rasters).
Therefore, we will transform the bathymetry data (depth_raster) and the economic exclusion zone data (eez) to match the coordinate reference system (crs = EPSG:9122) of sea surface temperature data (sst_rasters).
#define the target CRS (EPSG:9122) to be used in the transformations
target_crs <- crs(sst_rasters)
#transform the depth_raster to the target crs
depth_raster <- terra::project(depth_raster, target_crs)
#transform the eez to the target crs
eez <- st_transform(eez, crs = target_crs)
#create a summary list with all the files
spatial_file_list <- list(sst_rasters = sst_rasters,
depth_raster = depth_raster,
eez = eez)Now verify that all of the spatial files have the same coordinate reference system with the check_crs function before continuing with analysis:
#check that all the files have the same crs (target_crs)
all_crs_match <- sapply(spatial_file_list, check_crs, target_crs = target_crs)
#verify that all the files have the same crs
if (all(all_crs_match)) {
#if all the files have the same crs (all TRUE) then print
print("All spatial files have the target CRS")
} else {
#if all the files do NOT have the same crs (at least one FALSE) then print
print("The following files do not match the target CRS:")
print(names(spatial_file_list)[!all_crs_match])
}[1] "All spatial files have the target CRS"Since all the spatial files have the same coordinate reference system, we can move to the analysis process.
This section identifies locations that would be suitable for only bivalve farming:
To start, we will ensure that the sea surface temperature (SST) and the bathymetry data are compatible to process together
First, we need to calculate the average sea surface temperature by taking the mean across all of the rasters (from the years 2008 to 2012), and then convert this temperature from Kelvin to Celsius.
#find the mean sea surface temperature and store it as a single layer raster
sst_mean <- terra::mean(sst_rasters, na.rm = TRUE)
#now convert the temperature to Celcius by subtracting 273.15 to convert from Kelvin
sst_mean_celsius <- sst_mean - 273.15
#rename the layer for the celcisus raster
names(sst_mean_celsius) <- "mean_sst_c"Now we need to verify that the rasters are compatible to be combined. We will start by identifying the geometric properties of with the extract_spatial_properties function:
#create a summary list with all the files
sst_bath_file_list <- list(sst_mean_celsius = sst_mean_celsius,
depth_raster = depth_raster)
#identify the spatial properties of all the files in the list
spatial_properties_sst_bath <- extract_spatial_properties(sst_bath_file_list)
#generate a kable from the spatial properties extracted above
spatial_properties_kable_sst_bath <- spatial_properties_sst_bath %>%
kable("html",
col.names = c("Name", "Resolution", "Extent", "Origin", "CRS"),
caption = htmltools::tags$div(style = "text-align: center; font-size: 20px;",
htmltools::tags$strong("Geometric Properties Identification Results"))) %>%
kable_styling(full_width = FALSE, font_size = 14) %>%
row_spec(row = 0, bold = TRUE) %>%
kable_classic(html_font = "Times New Roman")
#print the kable
spatial_properties_kable_sst_bath| Name | Resolution | Extent | Origin | CRS |
|---|---|---|---|---|
| sst_mean_celsius | 0.04166, 0.04166 | ext(-131.98, -114.99, 29.99, 49.99) | -7.6e-06, -3.8e-06 | EPSG:9122 |
| depth_raster | 0.00417, 0.00417 | ext(-132, -114, 29, 50) | 0, 0 | EPSG:9122 |
It appears that the rasters have the same coordinate reference systems (since we transformed them to match above), but differ in their resolution, extent, and origin.
To ensure compatible geometric properties we will transform the data so both raster have the same resolution, extent, origin, and coordinate reference system. Since the bathymetry data has a larger extent than the sea surface temperature, we can crop this raster to match the extent and origin of the temperature raster.
#crop the bathymetry raster to match the extent/origin of the sst raster
depth_raster_crop <- terra::crop(depth_raster, sst_mean_celsius)Since the bathymetry data has a finer resolution than the sea surface temperature (SST) data, we will coarsen the bathymetry data by resampling the data using the nearest neighbor approach.
#resample the cropped bathymetry data using the nearest neighbor approach
depth_resampled <- terra::resample(depth_raster_crop, sst_mean_celsius, method = "near")Now we will verify that the transformations were successful and the resulting rasters have compatible geometric properties.
Now we need to verify that the rasters are compatible to be combined by ensuring the geometric properties have been succesfully transformed with the check_spatial_properties function:
#create a summary list with all the files
sst_bath_transform_file_list <- list(sst_mean_celsius = sst_mean_celsius,
depth_resampled = depth_resampled)
#identify the spatial properties of all the files in the list
spatial_properties_sst_bath_transform <- check_spatial_properties(sst_bath_transform_file_list)
#generate a kable from the spatial properties extracted above
spatial_properties_kable_sst_bath_transform <- spatial_properties_sst_bath_transform %>%
kable("html",
col.names = c("Name", "Resolution", "Extent", "Origin", "CRS"),
caption = htmltools::tags$div(style = "text-align: center; font-size: 20px;",
htmltools::tags$strong("Geometric Properties Verification Results"))) %>%
kable_styling(full_width = FALSE, font_size = 14) %>%
row_spec(row = 0, bold = TRUE) %>%
kable_classic(html_font = "Times New Roman")
#print the kable
spatial_properties_kable_sst_bath_transform| Name | Resolution | Extent | Origin | CRS |
|---|---|---|---|---|
| sst_mean_celsius | MATCH | MATCH | MATCH | MATCH |
| depth_resampled | MATCH | MATCH | MATCH | MATCH |
Since the transformations were successful, we can move forward with identifying suitable locations for oyster cultivation along the west coast.
Before we can identify potential finfish species that will be suitable options for marine aquaculture, we need to determine which locations are within the scope of the temperature and depth parameters for oysters.
To start, we will classify locations as “optimal” or “not optimal” for oyster cultivation based on the average sea surface temperature calculated above. This is done by creating a mask of the average sea surface temperature data (mean_sst_celcius) and assigning appropriate categorical levels using the suitable_locations function:
#define the minimum temperature
min_temp <- 11
#define the maximum temperature
max_temp <- 30
#generate a masked version of the temperature raster to identify optimal and not optimal locations
sst_mask <- suitable_locations(sst_mean_celsius, min_temp, max_temp)Now we classified locations as “optimal”, “not optimal”, or “land” for oyster cultivation based on the bathymetry data processed above. This is done by creating a mask of the bathymetry data (depth_resampled) and assigning appropriate categorical levels using the suitable_locations function:
#define the "minimum" depth (meters below sea level)
min_depth <- -70
#define the "maximum" depth (sea surface)
max_depth <- 0
#generate a masked version of the bathymetry raster to identify optimal and not optimal locations
depth_mask <- suitable_locations(depth_resampled, min_depth, max_depth)
#we will now rename the "not optimal" group 3 to be land (elevations above sea level)
levels(depth_mask) <- data.frame(id = c(1, 2, 3),
cats = c("not_optimal", "optimal", "land"))Once the physical condition rasters have been reclassified, we can combine the reclassified rasters (sst_mask and depth_mask) into a single raster that contains the optimal locations for oyster cultivation based on both parameters using the optimal_locations function:
#create a raster for the optimal geographic range for oyster cultivation based on the temperature and depth parameters
oyster_range <- optimal_locations(sst_mask, depth_mask)After generating a raster with the optimal location for oyster cultivation, we can overlay the economic exclusion zones (eez) data. Before combining this information, we must rasterize the eez data to make it compatible with the other SpatRaster objects. We will create an individual raster object for each region (Southern California, Central California, Northern California, Oregon, and Washington) of the economic exclusion zones as well as a raster stack which contains each region as a layer using the rasterize_geometries function:
#rasterize the regions and create a raster stack with all of the regions
eez_results <- rasterize_geometries(eez, oyster_range, "region")
#save the list of region rasters
eez_regions <- eez_results$rast_list
#save the raster stack
eez_stack <- eez_results$rast_stack
#save the list of boundary boxes for each region
eez_bbox_list <- eez_results$bbox_listNow we can visualize the suitable locations for oyster cultivation in each region of the economic exclusion zone along the west coast of the U.S. Let’s start by making a color palette to visualize each of the regions within the west coast EEZ.
#create a palette with 5 colors for the regions in the EEZ
region_palette <- brewer.pal(n = 5, name = "BuPu")
#assign the colors in the palette to region names
region_names <- c("Washington",
"Oregon",
"Northern California",
"Central California",
"Southern California")
#assign the colors in the palette to the list of names in region_names
region_palette_named <- setNames(region_palette, region_names)
#create a copy of the eez dataframe with a color that corresponds to each region in the region_palette_named
eez_color <- eez_caps %>%
mutate(color = region_palette_named[region])Now we can create an interactive map to visualize the optimal locations for oyster cultivation within each EEZ region:
#create a map that visualizes the suitable oyster cultivation locations and eez regions
tmap_mode("view")
#create an object to use with htmltools for better visualization
oyster_range_map <- tm_shape(oyster_range) +
tm_raster(title = "Suitability",
palette = c("not_optimal" = "lightgrey", "optimal" = "tomato3"),
labels = c("Not Optimal", "Optimal")) +
tm_shape(eez_color) +
tm_fill(col = "color",
popup.vars = c("EEZ Region" = "region"),
alpha = 0.5,
title = "EEZ Regions") +
tm_layout(main.title = "Suitable Oyster Cultivation Locations and EEZ Regions",
main.title.size = 1.5,
main.title.fontface = "bold",
main.title.fontfamily = "Times New Roman") +
tm_scale_bar(position = c("left", "bottom"))
#add a title for the interactive view map
map_title <- tags$h1("Suitable Oyster Cultivation Locations",
style = "text-align: center; font-family: Times New Roman; font-size: 24px; font-weight: bold;")
#combine the title and the map
oyster_range_map <- htmltools::tagList(map_title, tmap_leaflet(oyster_range_map))
#print the map
oyster_range_mapThis section identifies locations that would be suitable for BOTH bivalve and finfish aquaculture:
Start by loading in the csv with the metadata for the 26 potential finfish species identified in the background section:
#read in the csv for the potential species
potential_species <- read_csv(species_file)
#subset the data to only the species common name and scientific name
potential_species_tbl <- potential_species %>%
select(1:2) %>%
mutate(number = paste0(row_number())) %>% #add a number based on the row number in the df
select(3,1,2)
#select only the first 13 rows of the table
potential_species_tbl_1 <- potential_species_tbl %>%
slice(1:13)
#select only the last 13 rows of the table
potential_species_tbl_2 <- potential_species_tbl %>%
slice(14:26)
#combine the two tables
potential_species_comb <- cbind(potential_species_tbl_1, potential_species_tbl_2)
#tidy up the data
potential_species <- potential_species %>%
clean_names() %>%
mutate(common_name = str_replace_all(common_name, " ", "_") %>% str_to_lower(),
scientific_name = str_replace_all(scientific_name, " ", "_") %>% str_to_lower())Then print out a table with a list of the potential species being considered.
#print a kable of the potential species names
potential_species_kable <- potential_species_comb %>%
kable("html",
col.names = c("Number", "Common Name", "Latin Name", "Number", "Common Name", "Latin Name"),
caption = htmltools::tags$div(style = "text-align: center; font-size: 15px;",
htmltools::tags$strong("Potential Finfish Species"))) %>%
kable_styling(full_width = TRUE, font_size = 10) %>%
row_spec(row = 0, bold = TRUE, italic = FALSE) %>%
column_spec(column = c(3, 6), italic = TRUE) %>%
kable_classic(html_font = "Times New Roman")
#print the kable
potential_species_kable| Number | Common Name | Latin Name | Number | Common Name | Latin Name |
|---|---|---|---|---|---|
| 1 | Arrowtooth Flounder | Atheresthes stomias | 14 | Pacific Ocean Perch | Sebastes alutus |
| 2 | Bocaccio | Sebastes paucispinis | 15 | Pacific Sardine | Sardinops sagax |
| 3 | Canary Rockfish | Sebastes pinniger | 16 | Pacific Whiting | Merluccius productus |
| 4 | Chinook Salmon | Oncorhynchus tshawytscha | 17 | Petrale Sole | Eopsetta jordani |
| 5 | Chum Salmon | Oncorhynchus keta | 18 | Pink Salmon | Oncorhynchus gorbuscha |
| 6 | Coho Salmon | Oncorhynchus kisutch | 19 | Rex Sole | Glyptocephalus zachirus |
| 7 | Dover Sole | Microstomus pacificus | 20 | Pacific Rock Sole | Lepidopsetta billineta |
| 8 | English Sole | Parophrys vetulus | 21 | Northern Rock Sole | Lepidopsetta polyxystra |
| 9 | Flathead Sole | Hippoglossoides elassodon | 22 | Sablefish | Anoplopoma fimbria |
| 10 | Lingcod | Ophiodon elongatus | 23 | Shortspine Thornyhead | Sebastolobus alascanus |
| 11 | Northern Anchovy | Engraulis mordax | 24 | Sockeye Salmon | Oncorhynchus nerka |
| 12 | Pacific Cod | Gadus macrocephalus | 25 | Widow Rockfish | Sebastes entomelas |
| 13 | Pacific Halibut | Hippoglossus stenolepis | 26 | Yellowtail Rockfish | Sebastes flavidus |
Now we will use the functions defined above to evaluate the compatibility of each potential finfish species in each EEZ region. For each species, we will generate a SpatRaster for each EEZ region whose area has been classified as “optimal” or “not optimal” for the cultivation of that species and oysters. We will then stack these SpatRasters to generate a SpatRaster stack for each species containing a layer for each region. Additionally we will generate a tibble which ranks each EEZ region for each potential finfish species based on the percent of optimal area in the region. Finally, we will also create a map for each region layer generated for each potential species. The functions are used in the below order:
Output: This workflow will generate three outputs:
#create an empty list to hold the output of the for loop
species_rast_list <- list()
#create a dataframe column with the scientific names from the potential_species_tbl
map_name_df <- potential_species_tbl[, 3]
#convert the map_name_df to a list
map_name_list <- as.list(map_name_df)
#bind the map_name_list on to the potential_species dataframe
potential_species_names <- cbind(potential_species, map_name_df)
#rename the new column to be tidyverse friendly
potential_species_names <- potential_species_names %>%
rename(species_name_caps = 9)
#create a tibble to hold the output of the for loop
region_rank <- tibble(species = character(),
region = character(),
optimal_area = numeric(),
total_area = numeric(),
percent_cover = numeric(),
rank = numeric())
#create a list to hold the species maps
species_map <- list()
#for every species in the potential species dataframe
for(i in 1:nrow(potential_species_names)) {
#extract the properly formatted species name from the map_name_list
species_name_caps <- potential_species_names$species_name_caps[i]
#save the species name as the scientific name
species_name <- potential_species_names$scientific_name[i]
#extract the minimum temperature for this species
min_temp <- potential_species_names$min_temp[i]
#extract the maximum temperature for this species
max_temp <- potential_species_names$max_temp[i]
#since the thermal range of tolerance was predicted from a model for some species, we have classified this in other columns (preferred_min_temp and preferred_max_temp) and left the min_temp and max_temp column empty
##so if the min_temp is NA we will use the preferred_min_temp instead
min_temp <- ifelse(is.na(min_temp), #if the min_temp is NA
potential_species_names$preferred_min_temp[i], #use the preferred_min_temp
min_temp) #if is.na(min_temp) = FALSE then use the min_temp already extracted
##or if the max_temp is NA we will use the preferred_max_temp instead
max_temp <- ifelse(is.na(max_temp), #if the max_temp is NA
potential_species_names$preferred_max_temp[i], #use the preferred_max_temp
max_temp) #if is.na(max_temp) = FALSE then use the max_temp already extracted
#the min_depth and max_depth in the dataframe are recorded as negative numbers to represent meters below sea level to align with the format in the bathymetry data which has a vertical datum of mean sea level
##however the function wants the min_parameter value which would correspond to the max_depth, so extract the maximum depth for this species
min_depth <- potential_species_names$max_depth[i]
##the function wants the max_parameter value which would correspond to the min_depth, so extract the minimum depth for this species
max_depth <- potential_species_names$min_depth[i]
#FUNCTION ONE A: generate a masked version of the sea surface temperature data using the suitable_locations function
temp_rast <- suitable_locations(sst_mean_celsius, min_temp, max_temp)
#FUNCTION ONE B:generate a masked version of the bathymetry using the suitable_locations function
depth_rast <- suitable_locations(depth_resampled, min_depth, max_depth)
#FUNCTION TWO A: determine the optimal range for this species using the optimal_locations function
optimal_rast <- optimal_locations(temp_rast, depth_rast)
#FUNCTION TWO B: determine the optimal range for this species AND the oysters using the optimal_locations function
potential_locations <- optimal_locations(oyster_range, optimal_rast)
#FUNCTION THREE: generate a raster stack with each layer in the stack being a region in the economic exclusion zone
results_list <- extract_regions(potential_locations, eez_stack, map_name_list)
#FUNCTION FOUR: extract the species raster stack from the results list
species_stack <- results_list$rast_stack
#rank the regions for each species
species_rank <- rank_regions(species_stack)
#add the species name to the tibble
species_rank <- species_rank %>%
mutate(species = species_name)
#add the new species_rank tibble to the full tibble
region_rank <- bind_rows(region_rank, species_rank)
#extract the list of region maps from results list
map <- results_list[["map"]]
#finish formatting the map
map <- map +
tm_shape(eez_color) +
tm_fill(col = "color",
popup.vars = c("EEZ Region" = "region"),
alpha = 0.5,
title = "EEZ Regions") +
tm_layout(main.title = expression(italic(species_name_caps)),
main.title.position = "center",
main.title.size = 1.5,
main.title.fontfamily = "Times New Roman") +
tm_graticules()
#add a title for the interactive view map
map_title <- tags$h1(species_name_caps,
style = "text-align: center; font-family: Times New Roman; font-size: 24px; font-weight: bold; font-style: italic;")
#combine the title and the map
map <- htmltools::tagList(map_title, tmap_leaflet(map))
#add the map to the species_map list
species_map[[species_name]] <- map
#save the species results list to the species rast list
species_rast_list[[species_name]] <- results_list
}Now let’s visualize some of the results from our analysis.
This section summarizes the results of the above analysis through a variety of visualizations:
First, we can generate a heatmap to determine which regions are most suitable for aquaculture farms for both oysters and the 26 potential finfish species.
#make sure the regions are in order from south to north
region_rank$region <- factor(region_rank$region,
levels = c("southern_california",
"central_california",
"northern_california",
"oregon",
"washington"))
#create a heatmap of the potential species and their ranking per region
ggplot(region_rank, aes(x = species, y = region, fill = percent_cover)) +
geom_tile(color = "black", size = 0.5) +
geom_text(aes(label = round(percent_cover, 1)), #show the percent cover to 1 decimal
color = "black", size = 2) +
scale_fill_gradient2(midpoint = 0, #make the midpoint 0 so we can make the 0 percent_cover values the "azure" color
mid = "azure",
high = "tomato3",
name = "Rank",
breaks = seq(0, 5, by = 1),
limits = c(0, 5)) +
scale_y_discrete(labels = c("washington" = "Washington",
"oregon" = "Oregon",
"northern_california" = "Northern California",
"central_california" = "Central California",
"southern_california" = "Southern California")) +
labs(title = "Species vs Region Rankings",
subtitle = "Heatmap of Percent Optimal Cover",
x = "Species",
y = "Region") +
theme(plot.title = element_text(family = "Times New Roman",
face = "bold",
size = 16,
hjust = 0.5),
plot.subtitle = element_text(family = "Times New Roman",
size = 14,
hjust = 0.5),
axis.title = element_text(family = "Times New Roman",
face = "bold",
size = 14),
axis.text.x = element_text(family = "Times New Roman",
size = 10,
angle = 45,
hjust = 1,
vjust = 1),
axis.text.y = element_text(family = "Times New Roman",
size = 10),
legend.title = element_text(family = "Times New Roman",
face = "bold",
size = 12),
legend.text = element_text(family = "Times New Roman",
size = 10),
legend.position = "right",
legend.box.background = element_rect(color = "black", size = 1))
Then we can determine which regions will be the most optimal for creating farms for oyster and a finfish species:
#subset the region_rank dataframe to only include the regions that are the best for cultivation
top_regions <- region_rank %>%
group_by(species) %>%
filter(rank == 1) %>%
ungroup()
#create a table with the top regions for each of the potential species
##excludes regions that have no number 1 regions
top_regions_kable <- top_regions %>%
kable("html",
col.names = c("Species", "Region", "Optimal Area", "Total Area", "Percent Cover", "Rank"),
caption = htmltools::tags$div(style = "text-align: center; font-size: 15px;",
htmltools::tags$strong("Top Ranked Regions for Each Species"))) %>%
kable_styling(full_width = TRUE, font_size = 10) %>%
row_spec(row = 0, bold = TRUE, italic = FALSE) %>%
column_spec(column = 1, italic = TRUE) %>%
kable_classic(html_font = "Times New Roman")
#print the kable
top_regions_kable| Species | Region | Optimal Area | Total Area | Percent Cover | Rank |
|---|---|---|---|---|---|
| sebastes_pinniger | washington | 3224.738 | 67478.12 | 4.778938 | 1 |
| oncorhynchus_tshawytscha | washington | 3224.738 | 67478.12 | 4.778938 | 1 |
| oncorhynchus_keta | washington | 3224.738 | 67478.12 | 4.778938 | 1 |
| oncorhynchus_kisutch | washington | 3224.738 | 67478.12 | 4.778938 | 1 |
| hippoglossoides_elassodon | washington | 2842.345 | 67478.12 | 4.212248 | 1 |
| engraulis_mordax | washington | 3209.864 | 67478.12 | 4.756897 | 1 |
| sardinops_sagax | washington | 3224.738 | 67478.12 | 4.778938 | 1 |
| merluccius_productus | washington | 3180.263 | 67478.12 | 4.713028 | 1 |
| oncorhynchus_gorbuscha | washington | 3224.738 | 67478.12 | 4.778938 | 1 |
| sebastes_flavidus | washington | 2842.345 | 67478.12 | 4.212248 | 1 |
Finally, we can examine one of the regional maps generated for a potential finfish species. We will visualize the maps generated for the Sebastes pinniger (Canary Rockfish) as an example of the maps created for all of the potential finfish species:
#examine the maps generated for the sebastes_pinniger species
tmap_mode("view")
#extract the species map from the map list
sebastes_pinniger_map <- species_map[["sebastes_pinniger"]]
#print the map
sebastes_pinniger_mapThese results show that for all of the potential finfish species chosen, the Washington EEZ is the most suitable for paired oyster aquaculture. However, further parameters for each species should be examined more closely for potential areas with the Washington region before moving forward with development.